/*
 * @file
 * @date 18.05.2010
 * @author Ivan Babushkin (jogrms@gmail.com)
 */

#include <model.h>
#include <default.h>
#include <localtypes.h>
#include <stdlib.h>
#include <math.h>
#include <memory.h>
#include <stdio.h>
#include <omp.h>
#include <sys/time.h>
#include <time.h>

unsigned int arraySize;
int i, j;
float r, m, rx, ry, r2, r4, r8, r14, a, ax, ay, aFactor, sigma6m2;
Model* model;

int timeDiff(struct timeval *t1, struct timeval *t2) {
	int secDiff = (t2->tv_sec - t1->tv_sec);
	return secDiff * 1000000 + (t2->tv_usec - t1->tv_usec);
}

void* modelThreadRoutine(void* arg) {
	struct timeval lastTime, currentTime;
	int iteration;
	
	iteration = 0;
	gettimeofday(&lastTime, NULL);
	struct timespec waitDelay;
	waitDelay.tv_sec = 0;
	waitDelay.tv_nsec = 1000000000l;
	while (1) {
	while (model->running) {
		modelUpdate();
		iteration++;
		
		gettimeofday(&currentTime, NULL);
		if (timeDiff(&lastTime, &currentTime) > 1000000) {
			lastTime = currentTime;

			model->iterationRate = iteration;
			iteration = 0;
			printf("%d ips\n", model->iterationRate);
		}
	}
	nanosleep(&waitDelay, NULL);
	}
}

void readParams() {
	model->params.atomCount = ATOM_COUNT; // FIXME
	model->params.deltaT = DELTA_T;
	model->params.sigma = SIGMA;
	model->params.epsilon = EPSILON;
	model->params.atomMass = ATOM_MASS;
	model->params.boxSize = BOX_SIZE;
	model->params.frameDelay = FRAME_DELAY;
	model->params.stepFactor = STEP_FACTOR;
	model->params.tempFactor = TEMP_FACTOR;
	model->params.numThreads = NUM_THREADS;

	int e;
	int numi;
	float numf;
	char c;
	char* s = calloc(80, sizeof(char));
	e = scanf("%c", &c);
	while (e != EOF){
		if (c == '#') {
			while ((c != '\n') && (e != EOF)) {
				e = scanf("%c", &c);
			}
			e = scanf("%c", &c);
		}
		else {
			e = scanf("%s", &s[1]);
			s[0] = c;
			if (strcmp(s, "atom_count") == 0) {
				e = scanf("%d", &numi);
				if (e != EOF) {
					model->params.atomCount = numi;
				}
			}
			else if (strcmp(s, "frame_delay") == 0) {
				e = scanf("%d", &numi);
				if (e != EOF) {
					model->params.frameDelay = numi;
				}
			}
			else if (strcmp(s, "delta_t") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.deltaT = numf;
				}
			}
			else if (strcmp(s, "box_size") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.boxSize = numf;
				}
			}
			else if (strcmp(s, "epsilon") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.epsilon = numf;
				}
			}
			else if (strcmp(s, "atom_mass") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.atomMass = numf;
				}
			}
			else if (strcmp(s, "sigma") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.sigma = numf;
				}
			}			
			else if (strcmp(s, "step_factor") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.stepFactor = numf;
				}
			}			
			else if (strcmp(s, "temp_factor") == 0) {
				e = scanf("%e", &numf);
				if (e != EOF) {
					model->params.tempFactor = numf;
				}
			}		
			else if (strcmp(s, "num_threads") == 0) {
				e = scanf("%d", &numi);
				if (e != EOF) {
					model->params.numThreads = numi;
				}
			}
			else {
				while ((c != '\n') && (e != EOF)) {
					e = scanf("%c", &c);
				}
			}

			e = scanf("%c", &c);
		}
		while (((c == '\n') || (c == ' ')) && (e != EOF)) {
			e = scanf("%c", &c);
		}
	}
}

void initOMP() {
	omp_set_dynamic(0);
	omp_set_num_threads(model->params.numThreads);
}

Model* modelInit() {
	model = malloc(sizeof(Model));

	readParams();

	sigma6m2 = pow(model->params.sigma, 6) * 2.0f;
	arraySize = model->params.atomCount * sizeof(Point2f);
	m = model->params.atomMass * A_M_U;
	aFactor = 24.0f * model->params.epsilon *
		pow(model->params.sigma, 6) / m;

	model->p = malloc(arraySize);
	model->v = malloc(arraySize);
	model->a = malloc(arraySize);
	model->newa = malloc(arraySize);

	int side = ceil(sqrt((float) model->params.atomCount));
	float //step = model->params.boxSize / (float) (side + 1);
			step = model->params.sigma * model->params.stepFactor;

	memset(model->v, 0, arraySize);
	for (i = 0; i < model->params.atomCount; i++) {
		model->p[i].x = (float) (i % side) * step + step;
		model->p[i].y = (float) -((int) -((float) i / (float) side)) * step + step ;
		//model->v[i].x = 5.0e12f;
		//model->v[i].y = 5.0e12f;
	}

	initOMP();

	model->running = 0;

	return model;
}

void modelUpdate() {
	memset(model->a, 0, arraySize);
	
#pragma omp parallel shared(model, aFactor, sigma6m2) \
		private(i, j, rx, ry, r2, r4, r8, r14, a, ax, ay) 
	{

#pragma omp for
		
	for (i = 0; i < model->params.atomCount - 1; i++) {
		for (j = i + 1; j < model->params.atomCount; j++) {
			rx = model->p[j].x - model->p[i].x;
			ry = model->p[j].y - model->p[i].y;
			r2 = rx * rx + ry * ry;
			//r = sqrt(r2);
			r2 = 1.0f / r2;
			r4 = r2 * r2;
			r8 = r4 * r4;
			r14 = r8 * r4 * r2;
			a = aFactor * (r8 - r14 * sigma6m2);
			ax = a * rx;
			ay = a * ry;
			
//#pragma omp atomic
			model->a[i].x += ax;
//#pragma omp atomic
			model->a[i].y += ay;
//#pragma omp atomic
			model->a[j].x -= ax;
//#pragma omp atomic
			model->a[j].y -= ay;
		}
	}

#pragma omp for
	for (i = 0; i < model->params.atomCount; i++) {
		model->v[i].x += model->a[i].x * model->params.deltaT;
		model->v[i].y += model->a[i].y * model->params.deltaT;
		model->p[i].x += model->v[i].x * model->params.deltaT;
		model->p[i].y += model->v[i].y * model->params.deltaT;

		if (model->p[i].x < 0.0f) {
			model->p[i].x = -model->p[i].x;
			model->v[i].x = -model->v[i].x * model->params.tempFactor;
		}
		if (model->p[i].y < 0.0f) {
			model->p[i].y = -model->p[i].y;
			model->v[i].y = -model->v[i].y * model->params.tempFactor;
		}
		if (model->p[i].x > model->params.boxSize) {
			model->p[i].x = 2.0f * model->params.boxSize - model->p[i].x;
			model->v[i].x = -model->v[i].x * model->params.tempFactor;
		}
		if (model->p[i].y > model->params.boxSize) {
			model->p[i].y = 2.0f * model->params.boxSize - model->p[i].y;
			model->v[i].y = -model->v[i].y * model->params.tempFactor;
		}
	}
	}
}
